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Abstract 

A recently developed algorithm for nonlinear system per- 
formance analysis has been applied to an F16 aircraft to 
begin evaluating the suitability of the method for aerospace 
problems. The algorithm has a potential to be much more 
efficient than the current methods in performance analysis 
for aircraft. This paper is the initial step in evaluating this 
potential. 
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1 Introduction 

Real world systems are necessarily nonlinear and the 
range in which these systems can operate safely and 
within specifications is of paramount interest to both 
system designers and their eventual users. Typi- 
cally, controllers are designed and analyzed using lin- 
ear methods at different operating conditions. In the 
final design stages, the nonlinear system is analyzed 
through repetitive simulation to determine its actual 
performance. The currently accepted practice in in- 
dustry and elsewhere is to select a number of poten- 
tial operating parameters, select their operating range, 
and then randomly simulate the system [1], The pro- 
cess is known as Monte Carlo simulation. In order to 
ensure that the full range of the system’s behavior is 
covered, a very large number of simulations must be 
performed. In a typical scenario of an aircraft landing, 
5 weight conditions, 5 center of gravity locations, and 

2 flap settings are considered. The number of Monte 
Carlo simulations associated with this scenario, in or- 
der to ensure sufficient exploration of noise and dis- 
turbance space, is 5000. This scenario constitutes a 
small part of the operational envelope of an airplane, 
although an important one. As the number of parame- 
ters allowed to vary increases, the problem of checking 
performance robustness of a realistic system with this 
method grows exponentially in size and complexity. 

Any methodology that would provide for a reduction 
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in the number of simulations necessary to ascertain ro- 
bust performance of an airplane would be very useful 
in both control system development and in assisting 
the certification process. In recent work, we have de- 
veloped a numerical tool for nonlinear robust perfor- 
mance analysis that has the potential to reduce the 
number of simulations necessary to assess aircraft per- 
formance and thus ameliorate the exponential growth 
of the problem [5]. Emphasis was placed on creat- 
ing a computationally sound tool, requiring only in- 
formation usually available on the process being an- 
alyzed. This analysis tool only requires a simulation 
program for the plant. The development in many re- 
spects parallels that of robustness analysis for linear 
systems. We would like to ascertain the suitability of 
the developed algorithm for solving problems encoun- 
tered in aerospace applications. As the first test of 
the methodology, the algorithm is applied to an F16 
fighter aircraft executing a maneuver to determine the 
worst deviation from the trajectory. 

For linear time invariant (LTI) systems with com- 
plex, structured uncertainty, analysis of robust per- 
formance can be reduced to searching for the solution 
of a set of algebraic equations which give bounds on 
the achievable performance. One is thus able to find 
computationally efficient solutions, such as the power 
algorithm for the fi lower bound, without doing an ex- 
plicit parameter search involving repeated simulation. 
This works because the system is linear and the per- 
formance and uncertainty descriptions are chosen so 
as to give computationally attractive solutions, even 
for large problems. 

Performance analysis for nonlinear systems is diffi- 
cult due to the wide variety of behavior and structures 
which can occur. The algorithm summarized in this 
paper addresses the lower bound on the worst case per- 
formance. We will consider the problem of robust tra- 
jectory tracking: given a nominal trajectory for an un- 
certain, noisy nonlinear system, a feedback controller 
which stabilizes the trajectory and a description of the 
desired performance, find a lower bound on the worst 
case performance. The numerical tool we employ is a 
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power algorithm to solve a specific two point boundary 
value problem that is similar in computational nature 
to the power algorithm for the lower bound of fi, and 
so does not rely on classical descent techniques. 

The rest of this paper proceeds as follows. In Section 

2 a brief summary of the algorithm is given. In Section 

3 the F16 airplane model is described, the trajectory 
tracking problem is set up, and the results of applying 
the algorithm are discussed. This is an initial applica- 
tion of the recently developed algorithm to an aircraft 
problem. The ability of the algorithm to handle tabu- 
lar data, a relatively high number of parameters, and 
the class of nonlinearities present in aircraft, all typical 
characterizations of an aircraft problem, was unknown 
until this application. The last section concludes the 
paper with a discussion of results, issues encountered 
in this application, and new work. 


2 Brief Summary of Robust 
Tracking Problem Method 

Many nonlinear analysis problems of engineering in- 
terest can be reduced to a problem of tracking a nom- 
inal trajectory. Be it an airplane taking a turn or an 
idling engine going through a sudden change in load, 
the designer has in mind an appropriate path, to be 
completed in a finite predetermined time, and builds 
his control system accordingly. Since the real system 
is not exactly the one used for the design, and since it 
is also subject to noise, the system will not follow the 
intended trajectory. However, the design can still be 
considered successful if it remains close enough to it in 
an appropriate norm. 

In this paper we consider a restricted version of this 
problem. The original presentation of the algorithm’s 
development and a more detailed discussion can be 
found in [5]. Our performance measure will be the 2- 
norm of the error signal (i.e. the difference between 
the nominal and the actual trajectory). If needed, 
the error signal can be weighted by a multiplicative 
time function. Noise signals will be bounded in the 
2-norm. Unmodeled dynamics will be norm bounded 
operators. The only information available on these op- 
erators is their induced 2-norm. We will not restrict 
the operators to be causal. The system equations will 
be allowed to depend on a set of real parameters vary- 
ing in closed intervals. The initial conditions for some 
or all of the state variables will also be allowed to vary 
in given closed intervals. 

To simplify the notation we will work in the follow- 
ing with a system having one uncertain parameter, one 
unmodeled norm bounded component and one noisy 
input. All the results presented, however, generalize 


naturally. 

Let u be the noise signal perturbing the system, and 
let y be the error signal, that is the difference between 
the nominal and the actual trajectories. Denote by v 
the output of the uncertain dynamical block, and d the 
real uncertain parameter. The equations describing 
the system will then be 

x = f(x,u,v,8,t) 

y = g{x,u,v,8,t) 

z = h(x,u,v,8,t) 

with the following constraints 

1*1 < 1 

IM|2 = 1 

IM|2 = N| 2 . (1) 



Figure 1: Schematic representation for the robust tra- 
jectory tracking problem. 


Figure 1 gives a schematic representation of these 
equations. In order to determine how close the nomi- 
nal trajectory is being tracked we need to compute 


max 
|«5|<1, |M| = 1, 


IMI- 


( 2 ) 


The preceeding problem is a nonlinear constrained op- 
timization problem. It is in general non-convex and 
some of the optimization variables live in an infinite 
dimensional space. An exact solution is thus out of 
the question: we have to settle for upper and lower 
bounds. Note that this is true even if the system is 
linear: the complexity of the problem does not come 
just from nonlinearities in the system but from the 
nonlinear nature of the optimization index and the 
constraints. In the following section we briefly sum- 
marize the algorithm for computing a lower bound, 
based on the search for locally worst case signals. 
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2.1 Necessary Conditions for Worst 
Case Signals 

Any evaluation of the function ||j/|| 2 , for valid values 
of the parameters and signals is a lower bound on (2). 
So a simple way of getting lower bounds is through 
repeated simulation of the system for different val- 
ues of the uncertain signals in the model. This is 
at present the state of the art of nonlinear analysis 
as applied in industry: good simulation models are 
developed and designs are tested through extensive 
simulation, usually selecting the uncertain signals at 
random. This approach is practical, since it requires 
information from the plant that is usually available, 
and often gives reasonable results. The algorithm ap- 
plied in this paper improves on this approach without 
sacrificing in simplicity or in the generality of the in- 
formation required. Instead of simulating at random 
points, the algorithm looks for points that are good 
candidates for being local maximums and this search 
is conducted through a “power-like” method. The ro- 
bust trajectory tracking problem can be written as a 
standard Euler-Lagrange optimization problem. For 
a detailed development the reader is referred to [5]. 
The robust trajectory tracking problem is equivalent 
to optimizing the performance index 

J = j\dt=\\\yf (3) 

for the system verifying the differential equation 

x = f(x,u,v,x$) 

1 * 

x„ = —u u 

U 2 

f / * * \ 

= -(z z — V v) 

is = 0 , 

( 4 ) 

where 

V = g{x,u,v,xs) 
z = h(x,u,v,x$) 

with given initial conditions 

'T(Io) — XQ: x u ( i t 0 ^j — 0, ^a(Io) — 0? Xs(to^j — S (5) 

and final conditions 

x u (tf) = ^, x A (t f ) = 0. (6) 

So a set of signals u, v, and a parameter d achieve the 
worst case value of the performance index J only if 



r s = -i r 6 = 1 

A = 0, or < and or < and 

( A s (ti) <0 ( A s (ti) > 0. 

, (13) 

Remarks: Equations (7) and (10) describe a linear 
time varying dynamical system whose inputs are the 
(scaled) outputs of the original system. We will refer 
to this system as the adjoint or co-system. 

Equations (12) can be interpreted as an alignment 
condition between the outputs of the adjoint system 
and the inputs to the original dynamical system. Thus, 
these equations describe two dynamical systems inter- 
connected in a feedback loop. 

Equation (13) states that at an optimum, either the 
derivative of the performance index with respect to the 
value of the parameter is zero, or it is negative and the 
parameter is at the lower end of the interval or it is 
positive and the parameter is at the higher end of the 
interval. 

If we consider both the equations for the system, the 
co-system, and the alignment conditions together, we 
have a two point boundary value problem, i.e. a set 
of differential equations with boundary conditions at 
two distinct time instants. 

Several methods for solving the general two point 
boundary value problem have been devised and thor- 
oughly studied. (See for example [3], [2]). However, 
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the standard methods are based on gradient descent. 
In what follows we present a method to solve this par- 
ticular instance of the two point boundary value prob- 
lem that avoids the problems of gradient descent meth- 
ods. The algorithm is a generalization of the power 
algorithm for the lower bound of fi. In fact when ap- 
plied to linear systems the proposed algorithm reduces 
to the standard power algorithm for fi as described by 
Young and Doyle [6]. 


with direction z — v'j with the circle centered 

at (|y) z and with radius ||i>||. 

When solving the differential equations with a nu- 
merical integrator, we will obtain values for all the sig- 
nals at a finite number of time instants. The number of 
operations necessary to the signal updating described 
grows linearly with the number of time instants. 

The following iterative algorithm searches for fixed 
points of $, by evaluating it repeatedly. 


2.2 A Power Algorithm 

For a trajectory that meets the necessary conditions 
for a critical point, the Euler Lagrange conditions can 
be naturally separated into (i) a dynamical system 
with initial conditions only; (ii) a dynamical system 
with final conditions only; (iii) two sets of aligning 
conditions between the inputs and outputs of the two 
systems; and (iv) conditions relating the initial con- 
ditions of both systems. It is also important to note 
that the adjoint system depends on the trajectory of 
the original one. 

So, if the perturbations signals achieve the neces- 
sary conditions, the following composition of mappings 
yields the identity map: 

• Simulate the system along the current inputs. 


1. Simulate the system with « = 0, i’ = 0, x$ = 0. 
Use the time steps generated by the simulation 
routine as a time axis. 

2. Select random values for u along the time axis, 
(the signal is to be interpolated in between time- 
steps). Normalize u to fit the norm requirement. 
Set v° = 0, and x° s = 0. Set A^ = 1. 

3. Repeat 


(«*'+ V' +1 , z* +1 , AJ+ 1 , A* a +1 ) := $(«>•', 4 Aj,, A* a ) 


4. until 


i + 1 i + 1 i + 1 \i + 1 


A ’ A +1 )=(«>’' 


4al,a' a ) 


• Compute the co-system along the current trajec- 
tory, and simulate it backwards in time. 

• From the alignment conditions in (12) compute 
updated values for u, v, A u , and Aa- 

• Update the value of x$ with the following rule: 

X — T A sij'o') 

( -i x<-i 

x s = — !<X<! 

I 1 X > 1- 

Denote this composition by 

(u\ v\ S 1 , a; , Xi) = $(u 0 ,v 0 ,s°, A° , Aa). (14) 

Remarks: From the first equation in (12), using 
the old values of Aa and of the state trajectory we 
can compute uA u . Since A u is a scalar, and we know 
the norm of u we can separate this product into its 
components. 

From the second equation in (12) we can compute 
^(■|^) z — v'j Aa- Since we know the norm of v, and 

the value of z we can determine v and A a by 

intersecting the line passing through the origin and 


Remarks: If the algorithm converges, it converges 
to a fixed point of $ and thus to a set of signals that 
meet the necessary conditions for a critical point. 

In order to prove convergence we would have to 
prove that $ is a contraction around fixed points. That 
has not been proved even for the simpler case when the 
system is linear. Since we are interested in a prepon- 
derance of experimental evidence that this algorithm 
does in fact work with aerospace systems, we apply it 
to an F 16 fighter executing a climb and turn maneuver 
as our initial example. 

3 Application to an F16 

We want to determine whether the algorithm, summa- 
rized in the previous section, is suitable for aerospace 
applications. As a first step, the algorithm’s ability 
to handle a model that includes a number of nonlinear 
equations and tabular data with a relatively high num- 
ber of parameters that are allowed to vary, all charac- 
teristic of a typical aircraft, must be ascertained. 

The aircraft used in this example application is an 
F16. The aerodynamic model is a reduced version of 
the full model obtained in wind tunnel tests at NASA 
Langley in 1979 [4]. It consists of tabular data with 
typical interpolation routines and nonlinear equations 
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of motion. The engine model is that of an afterburning 
turbofan. The airplane model utilized in this applica- 
tion is defined for speed range of up to Mach 0.6 and 
angle of attack interval between -10 and 45 degrees. 
The model includes four traditional controls (elevator, 
aileron, rudder, and throttle) and thirteen states (ve- 
locity vector, attitude angles, angular velocities, nav- 
igational position, altitude, and engine power). Fur- 
thermore, the aerodynamic coefficients are built up in 
a traditional way and the equations of motion are full 
nonlinear flat earth equations. 

Several different trajectories were analyzed, but due 
to lack of space only one is presented as an illustrative 
example. The trajectory is a constant climb, constant 
g coordinated turn. The effective center of gravity (eg) 
location is set at 0.2c, the x-coordinate of the reference 
eg position normalized by the maximum aerodynamic 
chord length (c) of the wing, which makes the aircraft 
statically stable. The performance variables were cho- 
sen to be the turning radius error and the altitude 
error, both measured from a nominal condition char- 
acterized by absence of any exogenous disturbances as 
well as any uncertainty or parameter variation. The 
aircraft initiates the maneuver at 10,000 ft flying at 
500 fps. The F16 is then trimmed to climb at 50 fps 
while maintaining a 4.5 g coordinated turn. This is the 
aircraft’s nominal trajectory as illustrated in figure 2 
(solid line). 

During the maneuver the aircraft is subjected to at- 
mospheric turbulence in vertical, horizontal, and lat- 
eral directions modeled by Van Karman spectra and 
implemented by Dryden filters [1]. The magnitude of 
the allowable gusts is limited to 50 fps. In addition, 
seven parameters in the model are allowed to vary in- 
dividually on a closed interval. These parameters in- 
clude variation in eg position as well as uncertainty in 
the aerodynamic forces and moments along each axis. 
For the example presented in this paper the numerical 
values for the variations are as follows. Cg variation is 
on the interval between 0.195c and 0.205c. The aero- 
dynamic force coefficients are allowed to vary ±2.5% 
for Cx, ±1.5% for Cy, ±3% for Cz. The aerodynamic 
moment coefficients vary ±5% individually for rolling 
(Cl), pitching (Cm), and yawing (Cn) moments. 

The algorithm is asked to find the combination of 
parameters and wind gusts that produce the largest 
norm of the performance variable vector, i.e. turning 
radius and altitude error. The worst case combination 
produced by the algorithm gives the value of each of 
the parameters at the end point of the allowable in- 
terval of variation, eventhough the entire interval is 
searched. Numerically these are Cg at 0.195c, Cx at 
102.5% from nominal, Cy at 98.5%, Cz at 97%, Cl at 
95%, Cm at 95%, and Cn at 105%. The resulting 2- 


norm of the performance variables is 1260 ft, which is 
given for future comparison rather than physical mean- 
ing. 

The model simulation used by the algorithm was 
built up in a Simulink diagram, figure 6. The be- 
havior of the airplane under the worst case parameter 
variation selected by the algorithm is illustrated in fig- 
ures 2-5. The solid line in all the figures represents the 
nominal trajectory while the dashed line represents the 
perturbed trajectory. 

To compare with more traditional ways of evaluating 
nonlinear system behavior, Monte Carlo simulations 
were run. For each parameter the endpoints of the in- 
terval of variation were selected as allowable values. A 
system simulation with random turbulence subjected 
to the same restrictions as those imposed by the al- 
gorithm, i.e. random noise with normal distribution 
and limited energy input to Dryden filters, is run for 
each possible combination of parameter values, 128 in 
this case. For each of these parameter combinations 
10 simulations are performed. The resulting 2-norm 
of each simulation is plotted in figure 7. The figure 
shows the 2-norm of the performance vector for each 
of the simulations as well as the worst case 2-norm. A 
total of 1280 simulations were performed. 

As can be seen from figure 7, the 2-norm of the 
worst case parameter combination with atmospheric 
winds shaped by the algorithm is indeed larger than 
any combination of parameters with random atmo- 
spheric winds. The two vertical lines demarcate the 
interval that corresponds to the combination of pa- 
rameters selected by the algorithm as the worst case. 
Thus, the algorithm gives us a combination of param- 
eters that is particularly bad. While this combination 
is not unique, as is evident from the figure, it does pro- 
vide us with a better lower bound on the worst case 
behavior of the airplane for the allowable set of pa- 
rameter variations than the Monte Carlo method. In 
terms of computational efficiency, the worst case al- 
gorithm is at least four times faster than the Monte 
Carlo simulations in this particular case. It is impor- 
tant to note that both the algorithm and Monte Carlo 
simulations provide only a lower bound. It is possi- 
ble that even worse performance might be achieved 
by some other combination of parameters and atmo- 
spheric turbulence. There are currently no methods 
that compute the global maximum for a problem of 
this complexity. 

4 Conclusions and Future Work 

An application of a recently developed algorithm for 
robust performance assessment of nonlinear systems 
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to an F16 aircraft is presented in this paper. The 
algorithm successfully handles a 13 state nonlinear 
aero/propulsion model based on wind tunnel data in 
tabular form that is subjected to atmospheric turbu- 
lence and has a relatively high number of uncertain 
parameters. These results confirm the applicability of 
the nonlinear robust performance analysis method [5] 
to aircraft performance problems. 

In order to better match the needs of the aircraft 
designers and of the certification process, several is- 
sues must still be addressed. The first issue of interest 
is comparison between bounded energy and stochas- 
tic noise signals. The current practice in industry 
is to rely on stochastic performance measurements. 
Furthermore, for the certification the FAA requires 
stochastic performance measures. Future research will 
focus on incorporating the stochastics into the algo- 
rithm. 

Certification guidelines also establish probability 
distributions for uncertain parameters. We believe this 
can be incorporated into this method by subdividing 
the distribution curve into intervals and evaluating the 
worst case performance of the system while the pa- 
rameter falls into that interval is evaluated. We are 
currently studying this point. 

Although this algorithm cannot completely replace 
the Monte Carlo simulations necessary for the certifi- 
cation process, it does enhance the simulation results. 
It answers the question of how bad can performance 
of a system really get and under what circumstances 
would that behavior occur. This answer is extremely 
valuable during the development of control laws since 
it can be done cheaply in comparison to a large num- 
ber of Monte Carlo simulations and in parallel with 
linear robustness analysis. 
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Figure 2: Spacial view of the trajectory. Figure 4: Turning radius error 



-2000 -1500 -1000 -500 0 500 1000 1500 2000 

N-S Position, ft 


Figure 3: Ground track of the trajectory. 







Figure 6: Simulink diagram for the F16 robust performance analysis. 



Figure 7: Comparison of stochastic and worst case analysis. 
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